/** 
 * \file gdger.cu
 * \author Kyle E. Niemeyer
 * \date 10/06/2011
 *
 * Based on "dger.f" from BLAS.
 *
 */

////////////////////////////////////////////////////////////////////////

/** gdger performs rank 1 operation A = alpha * x * y' + A. BLAS level two.
 *
 * \param[in]   m       number of rows in matrix A
 * \param[in]   n       number of columns in matrix A
 * \param[in]   alpha   scalar alpha
 * \param[in]   x       vector, dimension (m-1)*|incx|+1
 * \param[in]   incx    increment between elements of x
 * \param[in]   y       array, dimension (n-1)*|incy|+1
 * \param[in]   incy    increment between elements of y
 * \param[in,out] A     matrix, dimension (lda, n)
 * \param[in]   lda     first dimension of A
 */
__device__ void gdger ( int m, int n, double alpha, const double *x, int incx, 
                        const double *y, int incy, double *A, int lda )
{
  
  // check input parameters
  if ( (m <= 0) || (n <= 0) || (lda < max(1,m)) || (incx == 0) || (incy == 0) || (alpha == 0.0) ) return;
  
  // start operations
  // elements of A are accessed sequentially with one pass through A
  
  if ( incy > 0 ) {
    uint jy = 1;
  } else {
    uint jy = 1 - (n - 1) * incy;
  }
  
  if ( incx == 1 ) {
    
    for ( uint j = 0; j < n; ++j ) {
      
      if ( y[jy] != 0.0 ) {
        double temp = alpha * y[jy];
        
        for ( uint i = 0; i < m; ++i ) {
          a[i + lda * j] += x[i] * temp;
        }
      }
      
      jy += incy;
    }
    
  } else {
    
    if ( incx > 0 ) {
      uint kx = 1;
    } else {
      uint kx = 1 - (m - 1) * incx;
    }
    
    for ( uint j = 0; j < n; ++j ) {
      
      if ( y[jy] != 0.0 ) {
        double temp = alpha * y[jy];
        uint ix = kx;
        
        for ( uint i = 0; i < m; ++i ) {
          a[i + lda * j] += x[ix] * temp;
          ix += incx;
        }
      }
      
      jy += incy;
    }
    
  }
  
} // end gdger